
setwd()



### read in data

amh <- read.csv("Amh_EPI.csv")
head(amh)



### crude protein yield models

a1 <- aov(amh$cp_yield_tha ~ amh$treatment)
summary(lm(a1))
anova(a1)
TukeyHSD(a1)

a1b <- aov(amh$cp_yield_tha ~ amh$treatment + amh$soilcol)
summary(lm(a1b))
anova(a1b)
TukeyHSD(a1b)
model.tables(a1b, "mean")

a1c <- aov(amh$cp_yield_tha ~ amh$treatment + amh$constraint)
summary(lm(a1c))
anova(a1c)
TukeyHSD(a1c)
model.tables(a1c, "mean")

a1d <- aov(amh$cp_yield_tha ~ amh$treatment + amh$topog)
summary(lm(a1d))
anova(a1d)
TukeyHSD(a1d)
model.tables(a1d, "mean")



### biomass yield models

a2 <- aov(amh$biomass_yield_tha ~ amh$treatment)
summary(lm(a2))
anova(a2)
TukeyHSD(a2)

a2b <- aov(amh$biomass_yield_tha ~ amh$treatment + amh$soilcol)
summary(lm(a2b))
anova(a2b)
TukeyHSD(a2b)
model.tables(a2b, "mean")

a2c <- aov(amh$biomass_yield_tha ~ amh$treatment + amh$constraint)
summary(lm(a2c))
anova(a2c)
TukeyHSD(a2c)
model.tables(a2c, "mean")

a2d <- aov(amh$biomass_yield_tha ~ amh$treatment + amh$topog)
summary(lm(a2d))
anova(a2d)
TukeyHSD(a2d)
model.tables(a2d, "mean")



### crude protein concentration model

a3 <- aov(amh$cp_pct ~ amh$treatment)
summary(lm(a3))
anova(a3)
TukeyHSD(a3)



### in vitro digestible concentration model

a4 <- aov(amh$ivd_pct ~ amh$treatment)
summary(lm(a4))
anova(a4)
TukeyHSD(a4)











